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ABSTRACT 


This  paper  uses  logic  programming  to  describe  numerical  integration  in  a general 
way.  Axioms  are  added  to  the  basic  logic  programs  to  define  specific,  known  numerical 
integration  algorithms.  Using  a slightly  different  formalization,  we  construct  logic 
programs  for  adaptive  Romberg  integration  and  for  a new  algorithm  for  adaptive  inte- 
gration. 

The  logic  programs  provide  a formal  basis  for  classifying  numerical  quadrature 
algorithms  (recall  Rice's  comment  that  there  are  1,000,000  useful  ones).  The  logic 
programs  are  also  more  understandable  than  the  corresponding  programs  in  conventional 
programming  languages. 

In  terms  of  logic  programming  there  are  several  novel  points  in  this  paper.  The 
data  type,  real  number,  is  not  defined  by  a recursive  constructor  function.  Termina- 
tion is  decided  by  error  bounds,  rather  than  by  reaching  some  basis  case  of  the  data 
constructor.  The  problem  itself  seems  to  demand  concurrent  execution  and  global 
variables  but  in  fact  is  solved  in  a linear  fashion. 


1.1  Review  of  Adaptive  Quadrature 


The  numerical  approximation  of  integrals  goes  back  to  the  early  Greeks  and  the 
approximation  of  tt.  An  excellent  survey  of  the  problem  and  modern  solution  methods 
is  found  in  Numerical  Integration  [2]. 

Ultimately,  all  numerical  approximations  take  the  form 


r 

J f (x)dx  = (b-a) £ wRf(xk) 


where  n Is  the  number  of  intervals,  and  where 

n-1 

£ "if1  • 
k*0 
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That  is,  the  approximation  is  a weighted  sum  of  function  values  over  the  interval 
of  integration.  The  problems  are: 

1.  choosing  the  sample  points  (x^) 

2.  choosing  the  weights  (wk) 

3.  knowing  how  accurate  the  result  is. 

Adaptive  algorithms  cluster  the  points  x^  within  regions  where  the  intergrand  is 
most  intractable.  The  strategies  for  doing  so,  and  the  resulting  programs,  have  tend- 
ed to  become  increasingly  opaque  and  complex  [4]. 

It  is  usual  to  present  only  deterministic  numerical  algorithms.  They  vary  by 
the  kind  of  approximation  used,  the  order  of  applying  approximations,  and  the  strategy 
for  termination.  If  we  broaden  our  view  to  include  nondeterministic  algorithms,  we 
can  arrange  them  hierarchically;  the  lower  in  the  hierarchy,  the  more  deterministic 
the  algorithm.  The  hierarchy  defines  families  of  algorithms  and  is  a natural  result 
of  using  logic  programming  for  specification. 

1.2  Review  of-  Logic  Programming 

Logic  programs  are  a subset  of  well-formed-formulas  (WFs)  of  first-order  predi- 
cate calculus  that  define  a function  or  relation  and  that  can  be  used  to  drive  the 
computation  of  a function  or  one  or  more  missing  values  of  a relation  [3,8]. 

The  WFs  of  logic  programs  are  restricted  to  the  clausal  form 

B «-  A]  a A2  ...  An 

where  the  A^ 's  are  predicates  of  the  form 

name ( parameter ,parameter2,. . • , parameter^ 

and  B is  either  a predicate  or  empty.  Each  parameter  is  either  a constant,  a variable, 
or  a simple  function  of  variables  or  constants.  The  variables  are  all  assumed  to  be 
universally  quantified.  It  has  been  shown  [1]  that  all  WFs  of  first-order  predicate 
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calculus  are  formal izable  within  this  restriction. 

The  procedural  interpretation  of  logic  programs  is  as  follows: 

«-  A,  A ...  A A 
1 n 

is  the  main  program  which  is  successfully  executed  when  each  is  known  to  be  true. 
Each  Ai  is  interpreted  as  a procedure  call. 

B <-  A-|  a A2  ...  An 

is  a procedure  definition.  B is  the  head  including  the  procedure  name  and  its  parame- 
ters; the  A^'s  are  the  procedure  body  and  correspond  to  procedure  calls  as  in  the  main 
program. 

A call  is  accomplished  by  matching  the  parameters  of  the  call,  A..,  to  a procedure 
head  B of  the  same  name,  binding  the  parameters  of  the  call  and  the  procedure  head. 
Each  call  has  the  effect  of  establishing  the  A. 's  as  subgoals. 

B - 

is  a basis  case  where,  having  been  called,  it  establishes  no  further  subgoals.  In  the 
predicate  sense,  B +■  is  an  assertion  that  8 is  true  for  the  parameters  given. 

A call  P(x-| ,. . . ,X|^)  is  carried  out  by  finding  a procedure  definition  of  the  form 

P(yr....yk)  *■  a q2  a ...  a qn 

and  then  unifying  (x^,y^)  1 5 i < k.  This  has  the  effect  of  assigning  values  to  the 
local  variables  of  the  procedure,  but  may  also  cause  values  to  be  passed  to  the  call- 
ing procedure. 

Example.  Consider  the  logic  program  for  the  factorial  function:^ 

1)  FACT(0,1 ) •*- 

2)  FACT(n+l,z)  «-  FACT(n.x)  a TIMES(n+l ,x,z) 

3)  «-  FACT(2,ANS) 


FACT(x.y)  is  true  if  and  only  if  factorial (x)  = y. 
TIMES(x,y,z)  is  true  if  and  only  if  x>y  * z. 


' ' ' ' . » * t 
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Clauses  1)  and  2)  are  two  separate  procedures  that  together  deftne  the  factorial 
function.  Clause  1)  is  the  procedure  that  terminates  the  computation  by  assigning 
the  output  to  be  1 if  the  input  is  0.  Clause  2)  is  the  procedure  that  reduces  a 
factorial  computation  to  a simpler  factorial  and  a multiplication.  We  assume  here 
that  TIMES  has  been  appropriately  defined  by  other  axioms  or  is  considered  a primi- 
tive function.  The  form  of  the  parameters  will  determine  which  procedure  definition 
is  actually  used  for  each  call.  Clause  3)  is  the  initial  call  and  serves  the  role 
of  forcing  a computation  of  factorial (2)  and  leaving  the  answer  in  variable  ANS. 

Let's  follow  the  computation.  FACT(2,ANS)  is  a call  to  a procedure.  It  cannot 
call  clause  1)  because  2 and  0 cannot  be  made  to  unify  since  they  are  two  distinct 
constants.^  So  a call  to  clause  2)  is  made.  FACT(2,ANS)  must  be  made  to  match 
FACT(n+l,z).  This  can  happen  if  we  bind  n to  1 and  ANS  to  z.  The  procedure  with 
local  variables  assigned  looks  like: 

FACT(2,ANS)  FACT(l,x)  a TIMES(2,x,ANS) 
i.e.  ANS  = 2-factorial (1).  This  gives  us  two  new  calls: 

FACT ( 1 ,x)  and  TIMES(2,x,ANS) . 


Normally  data  types,  e.g.  the  natural  numbers  used  here,  are  defined  constructive- 
ly, as  Peano  did,  e.g. 

NATNUM(O) 

NATNUM(n)  - NATNUM(s(n) ) . 


Then  clause  2)  would  be: 

FACT(s(n),z)  «-  FACT(n,x)  a TIMES(s(n) ,y,z) . 

The  above  can  either  be  interpreted  as  a shortcut  notation  for  Peano's  axioms,  or  for 
efficiency,  we  can  take  the  natural  numbers  and  a few  functions  on  natural  numbers  as 
primitives. 


-- .. — B 
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Let  us  first  consider  only  the  calls  to  FACT.  FACT(l.x)  cannot  call  clause  1)  because 
constants  0 and  1 will  not  unify;  clause  2)  is  called  again.  The  matching  process 


produces : 


FACT ( 1 ,x)  <-  FACT(O.x')  a TIMES(1 ,x' ,x) 


which.  In  turn,  produces  two  new  calls.  Considering  again  only  the  call  on  FACT: 
FACT(0,x')  calls  clause  1).  Clause  2 will  be  excluded  if  we  have  properly  defined 
the  natural  numbers,  i.e.  0 f n+1 . After  the  binding,  the  procedure  produced  from 


clause  1 is: 


FACT (0,1 ) 


which  is  unchanged,  but  the  binding  binds  1 with  x'  in  the  calling  procedure.  Since 
clause  1)  has  no  body,  we  simply  return.  Now  TIMES(l,l,x)  is  able  to  be  called  and 
binds  1 with  x,  giving  a fully  defined  procedure 

FACT (1 ,1 ) - FACT(O.l)  a TIMES(1 ,1 ,1 ) . 

In  addition,  x is  now  bound  in  the  previous  stage  and  as  we  return  we  find 

FACT (2 »ANS ) <-  FACT(l.l)  A TIMES(2,1 ,ANS) . 

TIMES  again  is  called  to  yield  ANS  = 2,  and  finally  FACT(2,2)  as  desired. 

We  use  arithmetic  expressions  in  the  parameter  lists  for  ease  of  readability. 

It  is  possible  to  compute  such  expressions  using  predicates,  e.g. 


/in-a  % 


could  be  computed: 


SUB(m,a,x)  a SUB(b,a,y)  a DIV(x,y,z)  a TIMES(z,e,e1 ) 


where 


SUB(u,v,w)  is  true  iff  u-v  = w 
DIV(u,v,w)  is  true  iff  u/v  = w 
TIMES(u,v,w)  is  true  iff  z*e  = e-j  . 

The  shorthand  notation  is  only  syntactic  sugar  and  does  not  change  the  essential 


character  of  logic  programs. 

The  implementation  of  logic  programming  systems  raises  some  points  not  directly 
relevant  to  this  paper,  i.e.,  how,  in  general,  to  select  procedures,  what  order  to 
evaluate  subgoals,  how  to  achieve  efficiency,  and  so  on.  To  some  extent  our  logic 
programs  reflect  these  considerations,  where  we  use  the  logic  itself  to  force  selec- 
tion and  order.  The  reader  interested  in  this  aspect  should  read  Warren's  description 
of  an  operational  logic  program  compiler  interpreter  [9]. 

2.1  Formal  Specification  of  Adaptive  Quadrature 

The  specification  is  a set  of  axioms  which  predetermine  the  desired  result.  The 
first  is  a formal  identity  from  the  calculus. 


f(x)dx  = f(x)dx  + jf(x)dx 

J J J 


From  the  above  we  can  derive  the  following  theorem: 


in  — 

Jf(x)dx  - y-j  < e-j  a |Jf(x)dx  - y2|  £ 
a m 

b 

- |Jf(x)dx  - (yj+y2)|  £ 6^+62 


Proof:  j s j <xa  Jtj  s y =»  |s+t|  - x+y* 

Formula  (1)  gives  the  basic  step  of  the  adaptive  routine,  reducing  the  whole 
problem  to  two,  presumably  simpler,  subgoals.  The  subdivision  process  terminates 
when  it  is  determined  that  a numerical  approximation  can  be  used  instead.  There 
are  many  ways  to  make  the  determination,  each  depending  upon  a pair  of  formulas, 
named  here  quad  and  err. 


* « / '"-V**  * • 
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For  example,  given  the  basic  trapezoidal  approximation 
trap(f ,a,b)  = (f(a)+f(b))- (b-a)/2 

we  may  use 


quad(f,a,b)  = trap(f,a,^-)  + trap(f.^-.b) 


and 


err(f,a,b;  = |quad(f,a,b)  - trap(f,a,b)| 
We  make  the  assumption: 

b 

|Jf(x)dx  - quad(f,a,b)  £ err(f,a,b) 


(2) 


While  neither  formula  (1)  nor  formula  (2)  is  true  in  general,  they  are  two  of  the 
assumptions  upon  which  numerical  quadrature  routines  are  based  and  therefore  have 
the  status  of  axioms  here. 

So,  given  interval  (A,B) , function  F,  and  error  bound  EPS,  we  would  like  to 
find  a y such  that 


B 

| J F(x)dx  - y 
A 


£ EPS 


(3) 


We  accomplish  this  by  subdividing  the  interval  and  error  bound,  as  described  in  (V) 
below  unless  err(f,a,b)  is  less  than  or  equal  to  the  error  bound  allotted  for  that 
interval,  in  which  case  we  evaluate  using  quad  and  then  return  that  answer  to  aid  in 
building  up  the  integral  for  the  next  inclusive  interval. 

The  structure  of  the  problem  can  be  seen  more  easily  perhaps  if  it  is  written  in 
predicate  form.  Define  the  predicate: 


I (a ,b,f , e,y)  « f(x)dx  - y 


U' 


£ e 





■■i 


* « * v X* 
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Then  consider  the  following: 
Logic  Program  1 : 
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I(a,b,f,e1+e2,y1+y2)  <-  err(f.a.b)  > 

(1 ' ) 

a I (a  ,m,f  ,e-j  .y^ 

a I (m , b , f , &2 yj 2 ) 

I(a,b,f,e,quad(f,a,b))  «-  err(f,a,b)  £ e 

(2') 

•<-  I (A, B,F, EPS, answer) 

(3’) 

(!')  is  the  recursive  case;  (2')  is  the  termination  case;  ( 3 ' ) is  the  main  program, 
i.e.  it  is  a call  to  I for  inputs  A,  B,  F,  EPS,  and  output  answer. 

The  above  logic  program  describes,  mathematically,  a large  class  of  integration 
algorithms.  However,  the  methods  for  splitting  the  interval  and  the  error  bound  are 
unspecified  and  the  computation  will  be  highly  nondetermini  Stic.  Specifying  the 
splitting  methods  will  give  us  particular  algorithms,  and,  in  general  more  determin- 
ism. In  addition,  we  may  not  even  get  closer  to  a solution  at  eacn  step.  Suppose, 
for  example,  that  we  always  chose  a "midpoint",  m,  outside  the  interval  (a,b).  Then 
instead  of  reducing  the  interval  size  at  each  step,  we  may  actually  be  increasing  it. 
Nothing  in  the  logic  program,  as  it  stands,  contradicts  this  behavior.  So,  let  us 
consider  some  ways  to  increase  the  efficiency  of  this  algorithm  by  further  specifica- 
tion. 

Given  a clause  C and  a property  p,  we  can  tighten  or  restrict  C by  p by  replacing 
clause  C by  clause  C *■  p.  For  example,  to  require  that  "midpoint"  m be  strictly  be- 
tween points  a and  b,  change  clause  (1')  to: 

O')  «-  (a  < m < b) 

which  is  equivalent  to  adding  the  restriction  to  the  list  of  subgoals  of  (1‘): 
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I(a.b,f,e1+e2,y1+y2)  «-  err(f,a,b)  > e^+e2 

a I (a  ,m,f  ,£^  ,y-| ) 
a I (m,b , f ,e2 ,y2 ) 
a (a  < m < b) 

Consider  the  particular  algorithm  in  which  the  intervals  and  error  bounds  are 
both  halved. 

Logic  Program  2: 

I(a,b,f,e,y1+y2)  «-  err(f,a,b)  > e 
a I(a,^2^-,f ,|,yi ) 
a I(^-,b,f,|,y2) 
plus  (2' ) and  (3* ). 

If  we  want  to  modularize  the  program,  so  that  there  is  a procedure  to  decide  how 
to  split  the  interval  and  error  bound,  we  could  write  the  program: 

Logic  Program  3: 

I(a,b,f,e,y)  «-  err(f ,a,b)  > e 
a I(a,m,f,e1»y1) 
a I (m ,b ,f ,e2 ,y2 ) 
a MID(a,b,m) 
a ER_DIV(a,b,m,e,ei ,e2) 

together  with  (2')  and  ( 3 ' ) , and  procedures  defining  MID  and  ER_DIV. 

If  we  are  dividing  the  intervals  and  the  error  bounds  proportionately,  as  in  the 
earliest  adaptive  quadrature  routines  [5],  then  we  could  define  ER_DIV: 

ER_DIV(a,b,m,e,e.j  ,e2)  ■*- 

_ /b-nu 

A e2  ‘ W*e 


ima. 


ft  f I 


v XT 


% »* 


If  we  define 


TO 


MID(a,b,m)  *-  m s 

then  with  this  definition  of  ER_DIV,  Logic  Program  2 and  3 are  equivalent,  as  can 
easily  be  proven.  Note  that  the  equality  symbols  represent  equality  and  not  assign- 
ment. Procedure  ERJDIV  expresses  a relation  that  holds  among  its  arguments;  the 
value  computed  depends  upon  which  arguments  are  already  bound.  E.g.  if  you  called 
ERJDIV  with  values  for  a,  b,  e,  e-j , then  it  could  compute  m and  eg.  For  more  dis- 
cussion of  varying  the  orientation  of  functions,  see  [3,7]. 

3.  Global  Algorithms 

The  usual  measure  of  efficiency  is  the  number  of  times  the  function  is  evaluated 
in  the  quadrature  summation.  Assumptions  1,  2 & 3 define  a class  of  adaptive  quadra- 
ture algoritnms.  While  it  would  be  nice  to  choose  an  optimal  member  from  the  class, 
it  is  not  possible.  Rice  [6]  claims  that  there  are  at  least  1,000,000  algorithms, 
each  optimal  in  its  own  domain.  Furthermore,  by  the  Mean  Value  Theorem,  one  function 
evaluation  is  enough  for  continuous  integrands,  if  only  one  knew  how  to  choose  the 
single  sample  point  In  the  interval. 

There  are  other  measures  of  efficiency.  In  particular  the  amount  of  memory 
needed,  which  can  be  applied.  Measures  of  this  sort  are  implementation  dependent. 

Nevertheless,  one  can  pick  a strategy  which  Is  both  intuitively  attractive  and 
performs  well  experimentally.  Since  at  each  splitting  of  a sub-interval  the  estimate 
of  the  truncation  error  over  the  whole  sub-interval  is  known,  we  may  choose  to  split 
the  pending  sub-interval  having  largest  error  estimate.  This  strategy  Is  new  as  far 
as  known  to  the  authors  and  Is  presented  below. 

The  next  two  algorithms  are  considered  global  because  there  is  a sense  of  having 
to  know  about  many  subgoals  at  once. 

In  this  algorithm,  we  continue  to  use  the  notion  of  subdividing  intervals,  but 


instead  of  distributing  the  error  bound  evenly  over  the  interval,  we  would  like  to 
reserve  most  of  the  error  for  the  trouble  spots  by  taking  advantage  of  the  accuracy 
In  the  easy  spots.  We  think  of  the  interval  having  been  broken  into  sub-intervals 
already,  and  associated  with  each  sub-interval  is  an  approximation  of  the  integrand 
and  an  error  approximation.  If  the  sum  of  the  error  approximations  is  greater  than 
the  total  error  bound  allowed  (e)  across  the  interval,  then  we  divide  the  sub-interval 
having  the  largest  error  estimate.  We  redefine  I to  use  this  method,  and  we  call  I 
as  before.  Logic  Program  4 is  a formal  specification  for  this  algorithm. 

Logic  Program  4: 

I(a,b,f,e,area)  «- 

G(f ,e,{(a,b,err(f ,a,b)) },err(f ,a,b) .area) 

G(f,e,intset,errsum,area)  «-  errsum  £ e 
a zQUAD(intset,f .area) 

G(f,e,intset, errsum, area)  «-  errsum  > e 

a GREATEST (intset, (a, b,e)) 
a MID(a,b,m) 
a intset'  = intset 

- {(a.b.e)} 

U {(a,m,err(f,a,m)),(m,b,err(f,m,b))> 
a errsum1  = errsum-e  + err(f,a,m)  + err(f,m,b) 
a G(f,e, intset' .errsum' .area) 

The  semantics  of  the  predicates  used  in  Logic  Program  4 are  given  in 
Table  1.  Logic  program  definitions  of  zQUAD  and  GREATEST  are  given  below. 


■1 
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Compute  the  sum  of  quad  applied  to  all  elements  of  an  INTSET. 

ZQUAD(*,f,0)  «- 

ZQUAD(T,f,area)  - T = T'U{(a,b,e)} 

a T - {(a,b,e)}  = T*  § 

a zQUAD(T',f,area') 
a area  = (area1  + quad(f,a,b)) 

Find  a sub-interval  with  the  greatest  error  approximation. 

GREATEST ({ (a, b,e) },(a,b,e) ) 

GREATEST(T,(a,b,e))  <- 

T = T,U{(a1,b1,e1)} 
aT"  < (a-j  »b-|  »e1 ) } = T' 
a GREATEST(T',(a2,b2,e2)) 
a MAX(  (a-j»b-j,e-j),  (a2,b2,e2)  ,(a,b,e)) 

where 

MAX((a1,b1,e1),(a2,b2,e2),(a1,b1,e1))  <^>62 

MAX((a-j  ,b-|  ,e-| ) , (a2,b2»e2) , (a2,b25e2))  +-  e^  < e2 

The  predicate  G will  compute  exactly  the  same  values,  in  the  same  order,  as  one 
of  the  (nondeterministic)  paths  of  Logic  Programs  1 or  3.  If  MID(a,b,m)  «-  m = 
then  G will  compute  exactly  the  same  values  as  one  of  the  paths  of  Logic  Program  2. 


These  first  two  conditions  are  needed  to  create  a choice  function  for  sets. 
They  mean  slightly  different  things  and  both  are  required. 
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predicate 

semantics 

6(f ,e,intset,errsum,area) 

Value,  area,  is  the  sum  of  the  integrals 
over  the  intervals  of  1 ntset  under  function 

f,  within  error  e.  Value,  errsum,  is  the 
sum  of  the  error  components  of  the  elements 
of  i ntset.  Errsum  is  redundant,  since  we 
could  recompute  it  each  time,  but  it  is  use- 
ful both  for  efficiency  and  clarity. 

2QUAD( 1 ntset , f , area ) 

✓ >» 

^ quad(f,a,b)  = area 

(a,b,e)  € intset 

GREATEST(intset,(a,b,e)) 

(a,b,e)  is  an  element  of  intset  having  the 
highest  third  component,  l.e.  error  approxi- 
mation. 

MID(a,b,m) 

as  used  before,  given  interval  (a,b)  this 
procedure  selects  a midpoint  m for  subdivis- 
ion. 

Table  1 


Program  4 can  be  improved  by  reducing  the  number  of  evaluations  of  f and  quad. 
Every  time  err(f,a,b)  is  called,  quad(f,a,b)  must  be  evaluated,  which  forces  evalua- 
tion of  f(a),f((a+b)/2)  and  f(b).  The  quad  and  f values  are  computed  redundantly  and 
are  discarded.  The  next  program  saves  f and  quad  values  in  such  a way  that  for  any 
points  a and  b,  f(a)  and  quad(a,b)  are  computed  at  most  once.  The  saved  values  make 
the  resulting  formulas  a little  more  cumbersome,  but  no  more  complex. 
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Four  other  changes  have  been  made  to  Program  4: 


1)  Intervals  are  expressed  as  left  endpoint  and  width  rather  than  left 
and  right  endpoints. 

2)  Instead  of  keeping  a set  of  intervals  and  finding  the  one  with  the 
greatest  error  estimate,  we  keep  a list  of  intervals  in  decreasing 
order  of  error  estimate.  The  symbol  ® denotes  an  infix  "cons"  of  an 
element  to  a list. 


3)  In  order  to  prevent  redundant  computations  of  f in  deriving  the  value 


of  quad,  we  have  to  fix  the  quadrature  rule;  we  have  chosen  the 
trapezoidal  rule,  but  we  could  have  chosen  others. 


I 4)  zQUAD  is  redefined  to  use  the  information  that  has  been  saved  in 

inti ist. 


(Logic  Program  5 Cont'd): 


G(f,e,intlist,errsum,area)  *■  errsum  5 e 
a ZQUAD( inti ist.area) 

G(f .e.intlist, errsum, area)  «-  errsum  > e 
a inti i st  = (a,h,fa,fm,fb,q,e)  ® list 
a fm1  * f(a  + h/4) 
a ^ * (fa  + 2fm.|  + fm)*h/8 

A e1  * |q1  - ((fa  + fm)-h/4)| 

a MERGE ( 1 i st , (a ,h/ 2 , f a , fm^  .fm.q^.e^.list') 
a fntg  = f(a  + 3h/4) 
a qg  = (fm  + 2^  + fb) • h/8 

a e2  = |q2  - ((fm  + fb) -h/4) | 

a MERGE(list'  ,(a  + hy^.h^.fm.fmg.fb^^) ,1  ist") 
a errsum1  = errsum  - e + e^  + e2 
a G(f,e,list",errsum' .area) 

ZQUAD(() ,0)  <- 

ZQUAD((a,h,fa,fm,fb,q,e)  ® list  , q + area) 

+■  ZQUAD(1  ist.area) 

MERGE(O.i.i)  - 

MERGE ( (a^ . , fa^ . fm^ . fb^ »q^  .e^j ) ® list  . (a2»h2ffa2»fm2»fb2*q2»62) » 
(^2  »ii2 » f^2 » fm2  »fb2  »q2  »e2 ) ® tfa^  »fm^  »f b'j  »q^ »e^ ) ® list) 

- e2  > ei 

MERGE ( (a^ . . fa^ . fm^ »fb^ .q^ . e^ ) ® list  . (a2.h2,fa2.fm2.fb2fq2fe2)» 
(fl^  *^1  *^®l * ^1  *^bi »Qi »®i ) ® list  ) 

- e2  < ei 

a MERGE(list,(a2,h2,fa2,fm2,fb2,q2,e2)  , list') 
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4.  Romberg  Integration 


One  can  improve  the  accuracy  of  the  integral  over  a given  set  of  sample  points 
by  either  resorting  to  higher  order  rules  or  by  applying  the  Romberg  convergence 
formulas  to  the  trapezoidal  rule  values.  Romberg  is  preferable  because  it  is  numeri- 
cally stable  [2].  Romberg  cannot  be  applied  over  arbitrarily  spaced  sample  points; 
in  its  usual  form  the  interval  of  integration  must  be  uniformly  divided  into  2 sub- 
intervals for  some  exponent  k.  This  property  has  generally  precluded  its  use  in 
adaptive  routines. 

In  the  algorithms  previously  described,  there  may  be  sub-intervals  which  are  in 

i/ 

fact  subdivided  exactly  2 times  for  some  k.  Romberg  can  be  applied  to  them  giving 
improved  approximations.  Then  all  of  the  approximations  can  be  summed  to  give  the 
final  result,  which  is  almost  always  an  improvement  over  the  unaccelerated  values. 
There  is  no  new  information  on  the  error  estimate  so  it  remains  the  same.  The  adap- 
tive Romberg  algorithm  is  the  same  as  Logic  Program  5 in  the  subdivision  of  the 
intervals,  but  different  in  the  computation  of  the  integral.  It  is  intuitively 
attractive  to  identify  the  powers  of  2 arising  naturally  out  of  the  subdivision 
process  with  those  needed  by  Romberg,  particularly  where  recursion  is  used  since  the 
calling  routine  contains  some  values  of  use  in  the  Romberg  rule.  The  algorithm 
presented  here  does  not  make  that  identification  and  as  a result,  is  better  able  to 

k 

apply  the  Romberg  rule.  The  reason  is  that  2 contiguous  sub-intervals  may  break 
across  the  recursive  structure  but  are  available  in  the  lists  used  here. 


Logic  Program  6: 

Change  the  termination  case  (clause  2)  of  Logic  Program  5 to: 
G(f,e,intlist,errsum,area)  «-  errsum  s e 
a SORT (inti ist.intlist* ) 
a EVAL(0,(),0,intlist' .area) 


> ^ i - 


rrr 


f: 


and  add  the  following  procedures: 

1)  Procedure  S0RT(x,y)  takes  list  of  intervals  x and  sorts  it  based  on  the  first 
value  of  the  7-tuple  representing  each  interval  to  form  list  y,  i.e.  the  left 
endpoint  of  the  interval.  The  actual  program  for  sort  is  not  included  here. 
It  is  a simple  program,  and  operates  similarly  to  MERGE. 

2)  Procedure  EVAL(left, mid, intsize, valist, area)  takes  a sorted  list  (valist)  of 
intervals  and  determines  if  there  are  any  sequences  of  sub-intervals  (mid) 
to  which  Romberg  acceleration  can  be  applied.  Any  adjacent  intervals  of 
equal  width  (intsize)  are  eligible.^  EVAL  collects  in  "mid"  contiguous  se- 
quences of  equal  sized  intervals  in  modified  form,  then  sends  them  to  REDUCE 
for  evaluation.  The  value,  left,  is  the  sum  of  the  Romberg  approximations 
of  the  integrals  of  the  intervals  to  the  left  of  mid.  The  value,  area,  is 
the  sum  of  left,  and  the  Romberg  approximations  of  the  integrals  of  the  in- 
tervals in  mid  and  valist. 

EVAL(left,mid, intsize, () .left  + x) 

«-  REDUCE(() ,mid,x) 

EVAL(left,mid, intsize, (a, h, fa, fm,fb,qamb,e)  ® valist, area) 

«-  intsize  = h 
a qab  = (fa  + fb) • h/2 
a ql  ist  = qab  <»  qamb 

a EVAL(left,mid  ® (h, fa ,fb,qlist) .intsize, val ist.area) 


EVAL(left, mid, intsize, (a, h, fa, fm,fb, qamb, e)  ® valist, area) 
a intsize  f h 
a qab  = (fa  + f b) • h/2 
a ql ist  = qab  ® qamb 

a EVAL(left  + x, (h, fa, fb, qlist) ,h, valist, area) 
a REDUCE((),mid,x) 

| • 


s The  equality  test  on  real  numbers  makes  sense  here  because  if  the  numbers  being 
compared  are  mathematically  equal,  they  are  computed  exactly  the  same  way  (repeated 

Hivicinn  hi/  9)  h onra  will  alcn  he  numorira  1 1 m onual 


3)  Procedure  REDUCL(1  Ist-j  ,list2,area)  takes  the  trapezoidal  approximations  pro- 
vided by  EVAL  and  computes  the  values  needed  to  apply  Romberg  acceleration, 
that  is  to  say  the  zero-level  Romberg  values,  by  pairwise  combining  contigu- 
ous intervals  of  the  same  step  size  until  a single  interval  is  constructed. 
The  formula 

Tk(a,h)  + Tk(a+h,h)  = Tk+1(a,h) 

where  i (a,h)  is  the  trapezoidal  approximation  on  interval  (a,a+h)  with  2 +1 
sample  points  provides  the  basic  computation  of  REDUCE.  The  final  result  is 
a list  of  the  form 

T°(a,h)  ® T1 (a,h)  e T2(a,h)  ...  ® Tn(a,h)  . 

REDUCE  pairwise  combines  the  intervals  of  list2  and  places  the  result  in 
1 ist-, . When  li$t2  is  exhausted,  1 ist1  and  list2  are  exchanged  and  the  pro- 
cess repeated  until  the  list  has  been  reduced  to  a single  interval.  If  the 
starting  list  is  not  of  length  equal  to  a power  of  two,  there  will  be  left- 
over intervals  when  the  intervals  are  pairwise  combined.  These  are  passed 
to  ROMS  at  that  time,  thus  are  not  considered  for  further  Romberg  accelera- 
tion. 

REDUCE(() ,() ,0)  - 

REDUCE ( i 1 i st , ( h ,fa ,f b ,ql i st) ,area1 +area2 ) 

«-  R0MB(1  ,qlist,area1 ) 
a REDUCE((),ilist,area2) 

REDUCE (ilist.O, area)  «-  ilist  f () 
a REDUCE((),ilist,area) 

REDUCE(11ist,(h,fa,fm,qlist^ ) 9 (h,fm,fb,qlist2)  ® ilist', area) 

«-  qab  = (fa  + fb)*h/2 
a SUMLIST(qlist^ ,qlist2»qlist3) 
a qlist  * qab  ® qllst^ 

a REDUCE(ilist  ® (2-h,fa,fb,qlist),nist',area) 
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where  SUML  1ST (ql  1 sti  .qlistg.qlist^)  takes  two  lists,  ql  1 st-j  and  qlist2,  and 
pairwise  sums  their  corresponding  components  to  form  qllst^. 

SUMLIST((), (),())  - 

SUML  I ST  (q  <8  L,r  <8  L‘,(q+r)  ® x) 

SUMLIST(L,L',x) 

4)  Procedure  ROMB(k,qlist,area)  takes  the  list  provided  by  REDUCE,  which  is 
also  the  set  of  values 

R°(a,h)  <8  R°(a,h)  ® RgCa.h)  ...  <8  R°(a,h) 

L XL  _ 

where  R*(a,h)  is  the  k— acceleration  on  interval  (a,a+h)  using  2 +1 
function  values.  We  have 

Rk  = ^^n’1  “ Rpl] )/ (4^-1 ) for  k>0 

and  R^(a,h)  is  the  desired  result.  The  parameters  of  ROMB  are  k,  the  level 
of  acceleration,  qlist,  the  list  of  values  at  that  state,  and  area,  the 
ultimate  answer. 

R0MB(k,(),r,r)  +-  length(r)  = 1 
R0MB(k,L,r^  <8  r2  <8  R.area) 

- y = (4kr2-r])/(4k-l) 
a R0MB(k,L  <8  y,r2  <8  R.area) 

R0MB(k,L, r.area)  «-  L^()  a length(r)  = 1 
a ROMB(k+l,(),L,area) 

5.  Correctness 

The  logic  programs  given  are  verifiable  since  each  procedure  is  an  easily  proved 
theorem,  (in  fact  many  are  theorems  found  in  numerical  analysis  texts).  For  example. 
Logic  Program  1 consists  of  two  procedures  that  correspond  to  the  two  theorems: 


*■— *.  « ' . — ■ 


» — . — M .. . . 


20 


Thm. 


(err(f ,a,b)  s e)  ■+  |Jf(x)dx  - quad(f,a,b)|  s e 
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formula  (2)  and  transitivity 


Thm. 


(err(f ,a,b)  > e) 
m 

a |Jf(x)dx  - y-,  | s 
a 

b 

a |Jf(x)dx  - y2|  < e2 
m 

21- 

Theorem  of  section  2. 

So  long  as  the  compiler/ interpreter  preserves  the  logical  meaning,  we  are  assured  of 
correctly  integrating. 

To  have  total  correctness,  we  must  also  prove  that  the  algorithms  terminate.  If 
the  looping  or  recursion  of  a program  Is  governed  by  counting  down  a natural  number 

| 


b 

■*  |Jf(x)dx  - (y1  + y2)  | s 
a 
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to  zero,  or  decomposing  a constructed  object  until  an  atomic  or  primitive  object  is 
reached,  then  you  can  easily  prove  termination.  In  the  examples  here,  the  reduction 
of  an  error  term  to  a specified  bound  causes  the  algorithms  to  stop.  The  proofs  of 
termination,  therefore,  are  more  difficult  here.  Proof  of  termination  of  Logic 
Programs  2 and  3 can  be  achieved  by  proving: 

Given  interval  (A,B)  and  function  F,  there  exists  an  interval  size  d such  that 
for  all  a and  b such  that  (A  s a < b 5 B)  a (b-a)  < d 

err(F,a,b)  < e-(jpj). 


5.  Conclusions 

We  have  shown  several  numerical  integration  algorithms  and  their  specifications 
in  logic.  The  logic  programs  are  unusually  concise,  compared  with  most  specifications 
of  numerical  integration  algorithms.  Each  procedure  must  represent  a theorem  in  nu- 
merical analysis.  Partial  correctness  is  thus  assured,  by  definition.  Termination 
is  trickier,  both  in  deriving  a mathematical  guarantee  of  termination  and  termination 
when  dealing  with  truncated  reals,  as  on  a computer. 

We  have  extended  logic  programs  to  use  on  real  numbers  and  error-bound  termina- 
tion, and  except  for  proof  of  termination  and  deciding  how  to  reduce  a problem  into 
subproblems,  this  use  of  logic  programming  is  not  significantly  different  from  logic 
algorithms  on  constructively  defined  data  types. 
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